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We extract the bulk viscosity of hot quark-gluon matter in the presence of light 
quarks from the recent lattice data on the QCD equation of state. For that purpose 
we extend the sum rule analysis by including the contribution of light quarks. We 
also discuss the universal properties of bulk viscosity in the vicinity of a second 
order phase transition, as it might occur in the chiral limit of QCD at fixed strange 
quark mass and most likely does occur in two-flavor QCD. We point out that a 
chiral transition in the 0(4) universality class at zero baryon density as well as the 
transition at the chiral critical point which belongs to the Z{2) universality class both 
lead to the critical behavior of bulk viscosity. In particular, the latter universality 
class implies the divergence of the bulk viscosity, which may be used as a signature 
of the critical point. We discuss the physical picture behind the dramatic increase of 
bulk viscosity seen in our analysis, and devise possible experimental tests of related 
phenomena. 
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I. INTRODUCTION 

Recently, it was pointed out that bulk viscosity in SU (3) gluo dynamics rises dramatically 
in the vicinity of the deconfmement phase transition [l| . The analysis was based on an exact 
sum rule derived from low-energy theorems of broken scale invariance, the high-statistics 



lattice data on the equation of state [2], and an ansatz for the spectral density of the 
correlation function of the energy-momentum tensor. The rapid increase of bulk viscosity 
has been associated with the fast growth of the thermal expectation value of the trace of the 
energy-momentum tensor (9)t = S — 3P, where £ is the energy density and P is the pressure. 
Very recently, bulk viscosity in SU(3) gluodynamics has been measured on the lattice by 
;he direct analysis of the correlation function of the trace of the energy-momentum tensor 
3I (see also 0]). The numerical results of the study j^] agree very well with the previous 
analysis [1]. Since the methods used in Refs. , and the lattice observables used as an 
input, are quite different, this agreement indicates that the dramatic growth (by about three 
orders of magnitude!) of bulk viscosity in the vicinity of the deconfmement phase transition 
is indeed a prominent feature of SU (3) gluodynamics. 

Is this behavior specific for the first order phase transition in SU(3) gluodynamics? What 
is the "microscopic" dynamics responsible for the growth of bulk viscosity? Can this growth 
be classified in terms of the general theory of critical phenomena? Can it be used in the 
experimental studies to isolate the signatures of phase transitions in heavy ion collisions? 
In this paper we address these questions by generalizing the previous analysis to the case 
of QCD with two light quarks and a strange quark, which is of direct relevance to heavy 
ion experiments. We will use as an input the very recent high statistics lattice data on 



the equation of state of QCD wit 
Columbia-Bielefeld Collaboration 



1 almost physical quark masses from the RIKEN-BNL- 

The paper is organized as follows. In section [Til we review the definition of bulk viscosity 
and the formalism relating it to the correlation function of the trace of the energy-momentum 
tensor based on Kubo's formula. In section II I II we discuss the low-energy theorems based 
on QCD renormalization group which constrain the zero-momentum correlation functions 
involving the trace of the energy-momentum tensor. In this Section we generalize the sum 
rule of Ref. [l| to full QCD with 2 + 1 flavors. In Section IIVI we briefly review the rel- 
evant lattice thermodynamics results, introduce and motivate our ansatz for the spectral 
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density, and use the sum rule to extract the bulk viscosity near the critical temperature. 
We find that the behavior of bulk viscosity in full QCD is qualitatively similar to the case 
of SU (3) gluodynamics. We then proceed to the discussion of uncertainties associated with 
our method, and estimate the error bars which should be associated with our result. In 
Section |V] we classify the critical behavior responsible for the rapid growth of bulk viscosity 
near the phase transition. We argue that at zero baryon density, this behavior belongs to 
the 0(4) universality class. The bulk viscosity in this case does not diverge at T = T c , but 
has a cusp. On the other hand, in the vicinity of the chiral critical point at finite baryon 
density, the critical behavior is in the Z(2) universality class. This means that the bulk vis- 
cosity should diverge at T = T c . This behavior should manifest itself in heavy ion collisions 
through the decrease of average transverse momentum of produced particles, accompanied 
by the increase in total multiplicity. This is due to both the increase in entropy associated 
with a large bulk viscosity, and the associated quenching of the transverse hydrodynamical 
expansion of the system ( "radial flow" ) . Finally, we summarize our findings and discuss the 
qualitative physical picture which arises from our analysis. In particular, we argue that the 
rapid growth of bulk viscosity favors the scenario of "soft statistical hadronization" in which 
the expanding system hadronizes at the phase transition by producing a large number of 
color screening soft partons. 



II. KUBO'S FORMULA FOR BULK VISCOSITY 

According to Kubo's formula the shear 77 and the bulk ( viscosities are related to the 
correlation function of the stress tensor 9ij(x), i,j = 1,2,3, as 

t](uj) ^SuS km + S im 5 k i - -SikSirr^j + ((u)6 ik 6i m 

= - lim / d 3 x /°°dte i ( wt - kr 3([0 ifc (t,r),0, m (O)]> (1) 
w J Jo 

Contracting the i, k and /, m indices gives the bulk viscosity as a static limit of the correlation 
function of the trace of the stress tensor 

C = - lim - / dt d 3 re^ ([9 u (x), 0^(0)]) . (2) 

This formula can be written in terms of Lorentz-invariant operators if we recall that the en- 
semble average of the commutator of the Hamiltonian H with any operator D in equilibrium 
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is given by 

d 3 x 9qo(x),0 



([H,0)U = i(^ =0. (3) 



eq \ / eq 



Thus we can write (T5]) as 



C = - lim - / dt d 3 re^ ([^(x),^(0)]> . (4) 

In the following we are going to calculate the bulk viscosity using the low energy theorems 
involving the Green's functions of the trace of energy-momentum tensor 9(x). Therefore, it 
is convenient to re-write Kubo's formula (|1J) using the retarded Green's function as 

ii r°° f 11 i i 

C = - lim - / dt d 3 r e wt iG R (x) = - lim - iG R (uj, 0) = — lim - lmG R (uj, 0) . (5) 
9 uj J Q J 9 ^->o u 9 w->o ijj 

The last equation follows from the fact that due to P-invariance, the function lm.G R (u, 0) is 
odd in uj while ReG R (u>, 0) is even in to. Note that the non-vanishing ( implies the existence 
of a massless excitation in the spectral density p. 
Let us define the spectral density 

p(uj,p) = --lmG R (u;,p). (6) 

7T 

Using the Kramers- Kronig relation the retarded Green's function can be represented as 

^ 1 [<*> ImG R (u,j?) ^ _ r P(u,P) 



G K (u,p) = - K -^ldu= KV du (7) 

J-oo U — UJ — IE J _ OQ UJ — U + IE 

The retarded Green's function G R (uj,p) of a bosonic excitation is related to the Euclidean 
Green's function G E (uj,p) by analytic continuation 

G E {uj,p) = -G R {iuj,p), uj>0. (8) 

Using (j2J) and the fact that p(uj,p) = —p(—oj,p) we recover 

<8 = lim 0) = 2 f°° du . (9) 

In the next section we use the low energy theorems to relate (5 to the energy density and 
the pressure of hot hadronic matter. 
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III. LOW ENERGY THEOREMS AT FINITE T 



In a conformally invariant theory the trace of the energy-momentum tensor vanishes. 
Therefore, Eq. (TJ1) implies that the bulk viscosity is a measure for the violation of conformal 
invariance. QCD is a conformally invariant theory at the classical level. However, quantum 
fluctuations generate the scale anomaly which manifests itself in the non-conseravtion of the 



dilatational current s. 



Of + @g 



(10) 



where we adopted the shorthand notation F 2 = F afJ-u F^ u and q, q are three-vectors in the 
quark flavor space (i.e. summation over flavors is assumed). (3(g) is the QCD /3-function, 
which governs the behavior of the running coupling 



/'- 



(3(g) 



(ii) 



At the leading order (3 = —bg 3 /(4:ii) 2 . 

Although the scale symmetry of the QCD lagrangian is broken by quantum vacuum fluc- 
tuations, there is a remaining symmetry imposed by the requirement of the renormalization 
group invariance of observable quantities 8]. This symmetry manifests itself in a chain of low 
energy theorems (LET) for the correlation functions of the operator 8g{x) = ((3/2g)F 2 (x). 
These low-energy theorems entirely determine the dynamics of the effective low-energy the- 
ory. This effective theory has an elegant geometrical interpretation [9f; in particular, glu- 
odynamics can be represented as a classical theory formulated on a curved (conformally 
flat) space-time background 10J. At finite temperature, the breaking of scale invariance by 
quantum fluctuations results in 9 = £ — 3P ^ clearly observed on the lattice for £77(3) 
gluodynamics 2|; the presence of quarks |ll| including the physical case of two light and a 
strange quark 5, or considering large N c fljl does not change this conclusion. 

n n n n 

The LET of Ref. [8|, 19[ were generalized to the case of finite temperature in [15|, [16( 
(lattice formulation has been discussed very recently in Ref. [17J). The lowest in the chain 
of relations reads (at zero baryon chemical potential): 

9 A - (12) 



<i 4 x{T» G W,D(0))= [Tg^-d] (0)t 



where d is a canonical dimension of the operator O. In particular, 



f d'x(Te G (x),e G (o)} = (r-^-Aj 



(13a) 
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J d 4 x (T9 G {x), 9 F (0)) = (t-^ - 3) (9 F ) T , (13b) 

where Qq and dp are the contributions to the trace of energy- momentum tensor from the 
gauge field and the quarks correspondingly defined in (fit)]) . 
Using ([10D and ffT3aj) . ffT3bj) we obtain 

J d 4 x(T6(x),6(0)} = ( T ^- 4 ) (0o)t + 2(tj^-3\ (6 f )t + J d'x (T9 F (x), 6 F (0)} 

* { T WT ~ 4 ) {9)t + { T WT ~ 2 ) ( ^ f)t ' (14) 

where in the last line we neglected the (T8 F (x), F (0)) correlation function which is propor- 
tional to m 2 . We will see in the next section that the quark mass m can indeed be treated 
small parameter . 

To relate the thermal expectation value of (#g)t to the quantity (S — 3P)* computed on 
the lattice, we should keep in mind that 

{£ - 3P)* = (6) T - (9) , (15) 
i.e. the zero-temperature expectation value of the trace of the energy-momentum tensor 

(6>) = -4|e,| (16) 
has to be subtracted; it is related to the vacuum energy density e v < 0. Analogously, 

(0f)t = (mqq)* + (mqq) . (17) 

Since the lattice data we are going to use correspond to almost physical quark masses, 
we can use the PCAC relations to express the zero-temperature vacuum expectation value 
(mqq)o of the scalar quark operator through the pion and kaon masses M n , Mk and decay 
constants f n , fx- 

(mqq), = -M* fl - M 2 K f K (18) 
Eq. Qj implies that the l.h.s. of (fllj) is just the constant <3. Using (j!3ap . (jl3bl) and 



Note that the heavy quark terms decouple in the low-energy matrix elements of the energy-momentum 
tensor canceling against the heavy quark part of the beta- function in the Qq operator 
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(fI5]) ; (II7p . (IIB)) in the r.h.s. of (fT4|) we derive the following sum rule 

= Ts ( -i - 3 ) - 4(8 - 3P) + ( T^- - 2 ) (mqq)* + 16\e v \ - 6(mqq) (19) 



i „A „„s a 



T s ^ " 3 J - " 3P ) + " 2 J < m ^>* + 16 N + 6 ( M ' + M k ft) (20) 

where s is the entropy density and c s is the speed of sound. In deriving ( fl9l) we used basic 
thermodynamic relations to relate temperature derivatives of energy density and pressure 
to the entropy density, s = dP/dT, specific heat, c v = dS/dT and the velocity of sound, 
= dP/dS = s/cy. Eq. (1201) is the main result of our paper. 



IV. BULK VISCOSITY FROM THE LATTICE DATA 



A. Spectral density 

In order to extract the bulk viscosity ( from ( 1201 we need to make an ansatz for the 
spectral density p. At high frequency, the spectral density should be described by perturba- 
tion theory; however, the perturbative (divergent) contribution has been subtracted in the 
definition of the quantities on the r.h.s. of the sum rule (Tl9|) . and so we should not include 
the high frequency perturbative continuum p(u) ~ a^u A on the l.h.s. as well. Indeed, when 
the frequency u is much larger than the temperature of the system T, uj ^> T the spectral 
density p(uj) should be temperature- independent. Therefore the subtraction of the expec- 
tation value of the trace of the energy-momentum tensor performed on the lattice should 
remove the high frequency perturbative part of the spectral density. 

Note that the subtraction of the zero-temperature expectation value (9)o ( 116]) removes 
the ultra-violet (UV) singularity 1/a 4 (a is the lattice spacing) in (Tl5l) . so that the resulting 
quantity (£ — 3P)* is UV finite. Such a subtraction implies the removal of local contact 
terms in the renormalized correlation function of the trace of the energy-momentum tensor. 

In the small frequency region, we will assume the following ansatz 

p(^,0) _9C 
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which satisfies flS} and Substituting fl2T|) in (Tl9l we arrive at 

9a; C = T S Q^- 3 )-4(£-3P)+(tA_2^ (mqq)* + 16\e v \ + 6{M* f*+M 2 K f K ) (22) 

The parameter ujq = ujq(T) is a scale at which the perturbation theory becomes valid. 
On dimensional grounds, we expect it to be proportional to the temperature, ujn ~ T. We 



estimate it as the scale at which the lattice calculations of the running coupling [18| coincide 
with the perturbative expression at a given temperature (see discussion below). 

We note that the expression for bulk viscosity (1221) consists of a thermal part that can 
be determined through lattice calculations, and a vacuum contribution, which we fixed 
using PCAC relations for the quark condensates and using the gluon condensate value 
|e„| 1//4 = 250 MeV. With this the vacuum part is given numerically by 

16|e v |(l + - ■ 1.6) ~ (560 MeV) 4 ~ (3T C ) 4 . (23) 
8 

Using the lattice results for trace anomaly, {£ — 3P)/T 4 , the fermionic contribution to it 
as well as the square of the velocity of sound presented in Ref. [5j] we can determine the 
bulk viscosity. To be specific we use the data set obtained from simulations on lattices with 
temporal extent N T = 6 j^]. These are the lattice results closest to the continuum limit. 
Our results for the bulk viscosity in units of the entropy density are shown in Fig. [TJ We 
note that the contribution from the fermion sector of the trace anomaly is at least a factor 3 
smaller than the gluonic contribution. Once lattice calculations with physical quark masses 
are performed we expect this contribution to become even less important. In Fig. [ITright) 
we show results for (/s using 500 MeV < ujq < 1500 MeV. 



B. Uncertainties of the method 



Our analysis relies on lattice results for the equation of state. Recent calculations with 
0(a 2 ) improved actions yield quite consistent results [si, Q, however it should be noted that 
they are still not extrapolated to the continuum limit and the light quark masses used in the 
calculation are about a factor two larger than the physical values. Since the sum rule we use 
is exact (up to the terms quadratic in light quark masses), the main source of uncertainty is 
the ansatz for the spectral density. We have made the simplest possible assumption about 
its shape consistent with general physical requirements - in particular, the spectral function 
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gluonic: T>0 

fermionic: T>0 

vacuum -- -- 

total 
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(a) (b) 

FIG. 1: Bulk viscosity in the high temperature phase of QCD versus temperature in units of the 
transition temperature. Part (a) of the figure shows various contributions to 9luq(T)(^/Ts given 
in (|22p . In part (b) of the figure we show the bulk viscosity in units of the entropy density for 
ujq = 0.5, 1, 1.5 GeV (top to bottom) which reflects the uncertainty in the determination of this 
scale parameter. 



must be linear for small frequencies. We also expect that the spectral density entering our 
sum rule should vanish above a certain frequency u ^> T at which the spectral density 
becomes perturbative and temperature-independent. 

We have estimated this frequency using lattice results on the temperature dependence of 
the running coupling 18|, Il9j |; ujq is taken to be the inverse distance at which this running 
coupling becomes approximately temperature independent and runs according to the zero 
temperature /3-functions. Choosing this scale involves some uncertainty. Moreover, we have 
to admit that there is a sizable uncertainty in our ansatz. A different functional form of the 
spectral density would change the numerical value of extracted bulk viscosity. Fortunately, 
the temperature dependence of bulk viscosity is not sensitive to this uncertainty, and so our 
main conclusion about the rapid growth of this quantity near the phase transition is ro- 
bust. Moreover, we note that the spectral density extracted from the analyses of correlation 



functions on the lattice 



ice sl 



20] is quite similar to our ansatz. 
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model 


a 





7 


5 


0(4) 


-0.21 


0.38 


1.47 


4.82 


Z{2) 


0.11 


0.33 


1.24 


4.79 



TABLE I: Critical exponents of 3-d 0(4) [2j| using (3 and 5 as input and Z(2) [22| symmetric spin 
models. 



V. UNIVERSALITY 



A. The case of zero baryon density 

Let us discuss here a generic second order phase transition as it might occur in the chiral 
limit of QCD at fixed strange quark mass and most likely does occur in two-flavor QCD. 
The critical behavior of thermodynamic quantities is in general controlled by two external 
parameters, the reduced temperature and the baryon chemical potential. 

In the vicinity of the critical point the behavior of bulk thermodynamic quantities is 
governed by thermal (yt) and magnetic (yh) critical exponents, which characterize the scaling 
behavior of the singular part of the free energy density, 

f(t,h) = ~]nZ = b~ 1 f(V*t,V'*h) . (24) 

Here b is an arbitrary scale factor. It is expected that the chiral phase transition in QCD 
can be described by an effective, three dimensional theory for the chiral order parameter, 
which in the case of two-flavor QCD would amount to an 0(4) symmetric spin model. 

The scaling behavior of the specific heat is controlled by the critical exponents a = 
(2y t — l)/yt- Exponents relevant for a discussion of other susceptibilities as well as the quark 
mass dependence of these quantities are (5 — (1 — yh)/yt, 7 = yt/yh and 8 = 2/^/(1 ~ Vh)- 
Their numerical values for Z(2) and 0(4) symmetric spin models in three dimensions are 
given in Table [B 

Eq. f l24"l) can be used to extract scaling laws for various quantities, valid in the vicinity 
of the critical point. We are at present only interested in the scaling behavior of the specific 
heat as function of temperature for vanishing external field 

c v (t) = c t~ a + const , (25) 
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which is obtained from Eq. fl2M after taking two derivatives with respect to t and choosing 
b = t~ l l yt ] a constant c can be both positive and negative. In the 0(4) universality class the 
exponent a is negative. The specific heat thus will not diverge but will only have a cusp. 
From the relations given by Eq. ( |T9|) we conclude that similarly the bulk viscosity will not 
diverge but will have a maximum at T c and the velocity of sound will not vanish but will 
only attain a minimum. 



B. Chiral Critical Point 



The situation may be different at the chiral critical point, i.e. the second order phase 
transition point that might exist in the QCD phase diagram at non-vanishing chemical 



potential 



231 ] . If it exists, this critical point belongs to the universality class of the 3-d 



Ising model, which has a positive specific heat exponent a. The situation here, however, is 
a bit more complicated as the energy-like and magnetization-like directions of the effective 
Ising model do not coincide with the temperature and symmetry breaking (quark mass) 
directions of QCD, nor is the latter controlled by the baryon chemical potential. Derivatives 
of the partition function with respect to temperature, which give the specific heat, thus will 
usually be related to mixed derivatives with respect to the energy-like and magnetization- 
like direction of the effective Ising-Hamiltonian characterizing the universal behavior in 
the vicinity of the chiral critical point. When approaching this critical point on a generic 
trajectory typical in a heavy ion collision, say for fixed s/ub, the critical behavior thus will 
not be controlled by a but by j//3S which will give the dominant singular behavior 23j. This 
leads to an even more rapid divergence of the specific heat, c v ~ h" 1 ^ 5 with j/PS ~ 0.8 
and where h = A\T — T c \/T c + B\hb — A*c|/ f^c parametrizes the distance to the critical point. 

As should be clear from our discussion given in Section HVl this also implies that the bulk 
viscosity will diverge at the critical point and the velocity of sound will vanish. Numerical 
evaluation of bulk viscosity close to chiral critical point would require lattice data at finite 
Hb- The sum rule in this case would also get modified due to the addition of terms containing 
derivatives in the chemical potential. We do not attempt such a numerical analysis in the 
present paper. However, the arguments for the divergence of bulk viscosity in the vicinity 
of chiral critical point which we presented above are sufficiently general. 
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VI. SUMMARY AND DISCUSSION 

Our results for the bulk viscosity obtained by combining low energy theorems with lattice 
results for QCD with a physical strange quark mass and almost physical light quark masses, 
indicate that the behavior is qualitatively similar to the one observed previously for the case 
of SU(3) gluodynamics - the bulk viscosity rises dramatically in the vicinity of the phase 
transition. Therefore the increase of bulk viscosity near the phase transition is a sufficiently 
general phenomenon and is not specific only to the first-order phase transition in SU(3) 
gluodynamics. 

We have argued that at zero baryon number density and in the limit of vanishing light 
quark masses the bulk viscosity will develop universal critical behavior which is related to 
that of the specific heat and will belong to the 3-dimensional, 0(4) universality class if the 
transition is second order in the chiral limit. The bulk viscosity in this case does not diverge 
at (T,fi) = (T c ,fj, c ), but has a cusp. On the other hand, in the vicinity of the chiral critical 
point at finite baryon density, the critical behavior is in the Z(2) universality class. This 
means that the bulk viscosity should diverge at T = T c . In both cases the growth of bulk 
viscosity can be attributed to the excitation of a massless scalar mode responsible for long- 
distance correlations. The rather simple relation of the critical behavior of the bulk viscosity 
to the (static) universality classes of QCD phase transitions arises from the relation to bulk 
thermodynamic observables generated by the low energy theorems. This is in contrast to 
ihe dynamic critical universality classes discussed in the context of liquid- vapor transitions 
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Large bulk viscosity signals strong coupling between the dilatational modes of the system 
and its internal degrees of freedom. We note that while for most physical substances the 
bulk viscosity is smaller or at most of the same order as the shear viscosity, this is not always 
true. For example, 3 He in the vicinity of the critical liquid-vapor point exhibits the ratio of 
bulk-to-shear viscosities in excess of a million 24]. 



What are the implications of our results? Since the growth of bulk viscosity signals 
an increase in entropy, our results imply that the expansion of hot QCD matter close to 
the phase transition is accompanied by the production of a large number of soft partons. 
These partons are produced so that the expanding quark-gluon system can hadronize - the 
produced partons screen the color charges of the quarks and gluons originally present in the 
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system. We thus come to the picture of "soft statistical hadronization" - "soft" because the 
produced partons carry low momenta, and "statistical" because the hadronization pattern is 
unlikely to depend upon the phase space distributions of "pre-existing" partons - the memory 
about these distributions is largely erased by the produced entropy. One may speculate that 
the association of inherent entropy with the hadronization process is analogous to the "black 
hole hadronization" scenario, in which confinement is associated with an event horizon for 
colored partons 
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23| 



In the vicinity of the chiral critical point, the divergence of bulk viscosity should mani- 
fest itself in heavy ion collisions through the decrease of average transverse momentum of 
produced particles, accompanied by an increase in total multiplicity. This is due to both 
the increase in entropy associated with a large bulk viscosity, and the associated quenching 
of the transverse hydrodynamical expansion of the system ("radial flow"). 
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